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Solar p-mode oscillations are excited by the work of stochastic, non-adiabatic, pressure fluctu- 
ations on the compressive modes. We evaluate the expression for the radial mode excitation rate 
derived by Nordlund & Stein (2000) using numerical simulations of near surface solar convection. 
We first apply this expression to the three radial modes of the simulation and obtain good agree- 
ment between the predicted excitation rate and the actual mode damping rates as determined 
from their energies and the widths of their resolved spectral profiles. These radial simulation 



modes are essentially the same as the solar modes at the resonant frequencies, where the solar 
modes have a node at the depth of the bottom of the simulation domain. We then apply this ex- 



pression for the mode excitation rate to the solar modes and obtain excellent agreement with the 
low I damping rates determined from GOLF data. Excitation occurs close to the surface, mainly 
in the intergranular lanes and near the boundaries of granules (where turbulence and radiative 
cooling are large). The non-adiabatic pressure fluctuations near the surface are produced by small 
instantaneous local imbalances between the divergence of the radiative and convective fluxes near 
the solar surface. Below the surface, the non-adiabatic pressure fluctuations are produced pri- 
marily by turbulent pressure fluctuations (Reynolds stresses). The frequency dependence of the 
mode excitation is due to effects of the mode structure and the pressure fluctuation spectrum. 
Excitation is small at low frequencies due to mode properties - the mode compression decreases 
and the mode mass increases at low frequency. Excitation is small at high frequencies due to the 
pressure fluctuation spectrum - pressure fluctuations become small at high frequencies because 
they are due to convection which is a long time scale phenomena compared to the dominant 
p-mode periods. 

Subject headings: sun:oscillations- sun:p-modes- sun:convection- sun:numerical simulation 



Introduction 1948 

1975 

Two ideas for the source of p-mode excita- 1977 



Lighthill 1952; Stein 1967, 1968; Crighton 
Ando & Osaki 1977; Goldreich & Keeley 
Goldreich & Kumar 1990; Balmforth 1992; 



tion have been pursued - overstability (as in Musielak et al. 1994). We show here that it is the 

pulsating stars) and turbulent Reynolds stresses Pd y work of sto chastic, non-adiabatic pressure 

(as in jet noise) (Biermann 1948; Schwarzschild fluctuations which is the primary mode excitation 
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mechanism (Stein & Nordlund 1991; Bogdan et al. 
1993; Goldreich et al. 1994). Near the surface, the 
non-adiabatic gas pressure (i.e. entropy) fluctua- 
tions dominate. They are produced by radiative 
cooling at the solar surface which is not locally and 
instantaneously exactly balanced by convective 
heat deposition. Below the surface, non-adiabatic 
turbulent pressure (Reynolds stress) fluctuations 
dominate. They are produced by the turbulent 
convective motions. 

In a previous paper (Nordlund & Stein 2000) we 
derived an exact expression for the stochastic exci- 
tation rate of the radial solar p-modes by the PdV 
work of non-adiabatic gas and turbulent pressure 
fluctuations on the mode compression. We now 
use realistic numerical simulations of near surface 
solar convection (depth about 2.5 Mm) to evaluate 
this expression (Stein & Nordlund 1998). Because 
the largest entropy and pressure fluctuations oc- 
cur near the solar surface and because modes with 
frequencies in the 3-4 mHz range, where the ex- 
citation rate is largest, are confined near the solar 
surface, these near surface simulations include the 
primary excitation and damping processes. 

2. Mode Excitation: Formalism 

The rate of energy input to the modes can be 
calculated starting with the kinetic energy equa- 
tion for the modes (Nordlund & Stein 2000). Ne- 
glecting the viscous stresses, 
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After integrating this equation over depth, the flux 
divergence term only gives contributions from the 
end points and is negligible. The buoyancy term 
is small because mass is conserved so there is no 
net mass flux. The last term is the PdV work, 
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There are several contributions to this work. The 
displacement, £, has contributions from the modes 
as well as the random convective motions. The 
pressure, 5P, has coherent contributions from the 
modes and incoherent contributions from the ran- 
dom convective motions. Both coherent and in- 
coherent contributions can be further divided into 



adiabatic and non-adiabatic terms. The dominant 
driving comes from the interaction of the non- 
adiabatic, incoherent pressure fluctuations, 

(5P nad = (SlnP -T 1 S\np)P , (3) 

with the coherent mode displacement, 
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This is a stochastic process, so the pressure fluc- 
tuations occur with random phases with respect 
to the modes. Therefore one must average over all 
possible relative phases between them. The result- 
ing rate of energy input to the modes is (Nordlund 
& Stein 2000) 
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Here, SP^ is the time Fourier transform of the non- 
adiabatic total pressure. Av = l/(total time in- 
terval) is the frequency resolution with which SP^ 
is computed. £ w is the mode displacement eigen- 
function, which is typically chosen to be real for an 
adiabatic mode. In that case, taking the complex 
conjugate of the pressure is not necessary, but we 
retain it for generality. The mode energy is 



] J dr p £ (£f = MjfiiR) . (6) 



Here M u is the mode mass and V UJ (R) is the mode 
velocity amplitude at the surface. Eqn 5 is sim- 
ilar to the expression of Goldreich et al. (1994) 
(eqn. 26), but involves no approximations. Hav- 
ing the numerical simulation data, we can evalu- 
ate this expression exactly without having to make 
approximations in order to evaluate it analytically 



3. Simulations 

We simulate a small portion of the solar pho- 
tosphere and the upper layers of the convection 
zone, a region extending 6x6 Mm horizontally 
and from the temperature minimum at —0.5 Mm 
down to 2.59 Mm below the visible surface. We 



1 A detailed description of how eqn 5 is evaluated is given in 
the appendix. Its validity is tested by applying it to the 
three radial modes of the simulation domain (sec. 4). 
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solve the equations of mass, momentum and en- 
ergy conservation in the form: 

^ = -u-Vlnp-V-u, (7) 

du PI 

— = -u-Vu + g VlnP+-V-a(3) 

at p p — 

— = -u-Ve V-u 

at p 

+ Qrad + Q vise 1 

(9) 

where a is the viscous stress tensor, Q ra d is the 
radiative heating and Q V i SC is the viscous dissipa- 
tion. 

We use a non-staggered grid of either 125 2 cells 
horizontally x82 cells vertically or 63 3 cells. The 
spatial derivatives are calculated using third or- 
der splines and the time advance is a third or- 
der leapfrog scheme (Hyman 1979; Nordlund & 
Stein 1990). The code is stabilized by a hyper- 
viscosity which removes short wavelength noise 
without damping the longer wavelengths. 

A large fraction of the internal energy is in the 
form of ionization energy near the solar surface, 
so we use a realistic equation of state (including 
the effects of ionization and excitation of hydrogen 
and other abundant elements and the formation 
and ionization of H2 molecules). The pressure is 
found by interpolation in a table of P(ln p, e) and 
its derivatives which is calculated with the Upp- 
sala stellar atmosphere package (Gustafsson et al. 
1975). 

Radiative energy exchange is critical in de- 
termining the structure of the upper convection 
zone. Near the surface of the sun, the energy 
flow changes from almost exclusively convective 
below the surface to radiative above the surface. 
The interaction between convection and radiation 
near the solar surface determines what we observe 
and produces the entropy fluctuations that lead 
to the buoyancy work which drives the convection 
and that give rise to the non-adiabatic pressure 
fluctuations which excite the p-mode oscillations. 
Hence, the interaction between convection and ra- 
diation is crucial for both the diagnostics and the 
dynamics of convection. Since the top of the con- 
vection zone occurs near the level where the con- 
tinuum optical depth is one, neither the optically 
thin nor the diffusion approximations give reason- 



able results. We therefore include 3D, LTE, non- 
gray radiation transfer in our model. 

We simulate only a small region near the so- 
lar surface and must therefore impose boundary 
conditions inside the convectively unstable region. 
What happens outside our computational domain 
in principle influences the convective motions in- 
side. However, at the top boundary, the density 
is very low relative to the rest of the volume and 
hence whatever happens there has very little influ- 
ence on the convective motions. At the bottom, 
the incoming fluid is to a very good approximation 
isentropic and featureless, and hence carries little 
information. The unknown influence of external 
regions should therefore be small. This assertion 
is indeed confirmed by experiments with bound- 
aries located at different depths. 

The horizontal directions are taken to be peri- 
odic. In the vertical direction, we have a trans- 
mitting boundary at the temperature minimum 
(Nordlund & Stein 1990). This is achieved by 
a larger than normal zone at the top boundary. 
Across this zone we make the vertical derivative 
of the density hydrostatic, set the vertical deriva- 
tive of the velocity zero and hold the internal en- 
ergy at the top fiducial layer constant in time and 
space. At the bottom of the computational do- 
main, outgoing fluid goes out with whatever prop- 
erties it has. For incoming fluid, we adjust the 
pressure such that the net mass flux through the 
bottom boundary vanishes. (This ensures that 
there is no boundary work done on vertical oscilla- 
tion modes.) The pressure on the bottom bound- 
ary thus varies in time but is uniform over the hor- 
izontal plane. We damp fluctuations of the hori- 
zontal and vertical velocity of the incoming fluid, 
using a long time constant. Finally, we adjust the 
density and energy of the incoming fluid, at con- 
stant pressure, to fix its entropy (in both space 
and time). 

The ability to do a direct numerical simula- 
tion with the wide range of length scales matching 
the dimensionless parameters - Reynolds number, 
Raylcigh number and Prandtl number - of the so- 
lar convection zone is beyond the speed and mem- 
ory capabilities of current computers. Thus, our 
simulations are of the type called "Large Eddy 
Simulations". It is, however, possible to resolve 
the surface thermal boundary layer of the convec- 
tion zone and this we have done. Indeed, this is re- 
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quired to achieve results that agree quantitatively 
with solar observations (Stein & Nordlund 2000). 

The picture of convection that has emerged 
from the simulations is the following: Convection 
is driven by radiative cooling in the thin thermal 
boundary layer at the solar surface. It consists 
of cool, low entropy, filamentary, turbulent, down- 
drafts that plunge through a warm, entropy neu- 
tral, smooth, diverging, laminar upflow. Upflows 
must diverge as they ascend into lower density lay- 
ers in order to conserve mass. This divergence 
smooths out any variations in their properties that 
might arise. Only a small fraction of the fluid at 
depth reaches the surface to be cooled and form 
the cores of the downdrafts. Most fluid turns over 
within a scale height and is entrained by the down- 
drafts. These low entropy downflows are the site 
of most of the buoyancy work that drives the con- 
vection. (See Stein & Nordlund (1998) for more 
details.) 

We have made numerous comparisons between 
the predictions of the simulations and solar obser- 
vations: 

• The depth of the convection zone depends 
on the opacities that determine the temper- 
ature versus pressure (and hence entropy) 
stratification of the surface layers, the spec- 
tral line blocking, the convective efficiency of 
the superadiabatic layers immediately below 
the surface that determines the transition to 
the asymptotic adiabat, and the equation of 
state that determines the further run of tem- 
perature versus pressure through the convec- 
tion zone. Excellent agreement is obtained 
between the depth of the convection zone 
predicted by our numerical simulations and 
that inferred from helioseismology (Rosen- 
thal et al. 1999). 

• The cavity for high frequency modes is en- 
larged by turbulent pressure support and 
3D non-linear opacity effects which increase 
the average temperature required to pro- 
duce a given effective temperature in an in- 
homogeneous compared to a homogeneous 
atmosphere. The p-mode eigenfrequencies 
calculated from the mean simulation atmo- 
sphere are significantly closer to the ob- 
served mode frequencies than those for stan- 
dard spherically symmetric, mixing length 



models (Rosenthal et al. 1998, 1999). 

• The simulation granulation size spectrum 
and the distribution of emergent intensities, 
when smoothed by the point spread function 
appropriate for the Swedish Vacuum Solar 
Telescope on La Palma (which produces the 
best solar images available today), closely 
match the observations (Stein & Nordlund 
1998). 

• The width of photospheric iron lines (whose 
thermal speed is small) is a signature of 
the convective velocities. The net Doppler 
blue shift and asymmetry of spectral lines 
is a signature of the correlation between ve- 
locity and temperature fluctuations. Both 
these signatures agree closely with observa- 
tions (Asplund et al. 2000). 

This gives us confidence that the simulations are 
properly modeling the crucial properties of near 
surface solar convection. 

4. Mode Excitation: Simulation 

Three radial modes exist in our simulation and 
we first apply eqn (5) to these modes and com- 
pare the rate of work on the modes it predicts 
with the actual excitation rate dE^/dt = TE^ de- 
termined from the mode's energies and widths in 
the simulation. The modes can be clearly seen in 
the spectrum of the horizontally averaged, depth 
integrated kinetic energy (fig 1). Their properties 
are given in table 1. The lowest mode (with no 
zero crossings inside the computational domain) 
has a frequency of 2.57 mHz and a FWHM of 19 
/xHz. Hence a simulation significantly longer than 
14.5 hrs is required to resolve this mode. We use 
a simulation of 43 hours, with a spatial resolution 
of 63 3 . Snapshots were saved at 30 s intervals. 



Table 1: Mode Properties 
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IS.u 
(/iHz) 


E u 
(ergs) 


Q 


r 

(s- 1 ) 


dE/dt 
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3.4 x 10 26 
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2.6 x 10 2b 
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2.7 x 10~ 3 


6.1 x 10 21 



4 




Fig. 1. — Kinetic energy, horizontally averaged 
and integrated over depth. Three radial modes are 
clearly visible. Least squares, Lorcntzian, fits to 
the modes and linear fits to the background noise 
are superimposed. 

The eigcnfunctions of the velocity are calcu- 
lated by taking the time Fourier transform of the 
velocity. To get the real eigenmode, the trans- 
form of the velocity at the frequency of the mode 
is divided by its most common phase among all 
depths. To reduce the noise, we average the result 
over a frequency band, approximately equal to the 
FWHM of the mode, centered on the mode. The 
eigenfunctions are essentially the same as the so- 
lar model modes of Christcnscn-Dalsgaard using 
his spherically symmetric model S (Christensen- 
Dalsgaard et al. 1996) at the mode frequencies 
(fig 2), when the modes are normalized by the 
square root of the mode energy, or which is equiv- 
alent by their amplitude at the surface. Hence, we 
choose to use the solar model eigenfunctions, be- 
cause they are much denser (35 radial modes below 
the acoustic cutoff frequency instead of three) and 
because they are slightly smoother. Another way 
of looking at the modes is via their kinetic energy. 
Figure 3 is an image of the logarithm of the ki- 
netic energy as a function of depth and frequency. 
The three modes of the computational domain are 
clearly seen. Notice also that there is a continuous 
pattern of nodes and antinodes, with the nodes de- 
scending in depth as the frequency increases and 
new nodes starting at the surface. This pattern 
extends even into the propagating region above 
the acoustic cutoff frequency of about 5.3 mHz. 



Fig. 2. — Velocity eigenmodes of the simulation 
compared with those of a solar model (model S 
of Christensen-Dalsgaard et al. 1996). The modes 
are normalize by the square root of their energy 
(eqn 6). The solar modes that have nodes at the 
bottom of the simulation domain closely match the 
simulation modes. The simulation pi mode corre- 
sponds to the solar p 16 mode and the simulation 
Pi mode corresponds to the solar p 2 e mode. The 
simulation p$ mode lies above the acoustic cutoff 
frequency. 

The actual mode energy decay rate of the simu- 
lation modes, which on average is equal to their ex- 
citation rate, is obtained by multiplying the energy 
of each mode by its decay rate T = 27tA^fwhm, 
which is its radian line width and is obtained from 
the fit to the modes (fig 1). 

The total mode energy is the sum of the mode 
energies minus a fit to the background over the 
frequency range where the modes are significantly 
above the background multiplied by the area (36 
Mm 2 ) of the simulation domain (fig 1). 

Eqn (5) is used to predict the mode excita- 
tion rate. The non-adiabatic total (gas -I- turbu- 
lent) pressure fluctuations are calculated directly 
from the simulation by first averaging the gas pres- 
sure, turbulent pressure and density over horizon- 
tal planes at each saved snapshot. These are then 
interpolated to the Lagrangian frame at each time. 
The non-adiabatic pressure fluctuations are calcu- 
lated as in eqns Al and A2. The oscillation modes 
of the domain are essentially adiabatic and do not 
affect the non-adiabatic pressure fluctuations as 
can be seen from its featureless, noisy spectrum 
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Fig. 3. — Logarithm of kinetic energy as a func- 
tion of frequency and depth. The three resonant 
modes of the simulation stand out as maxima in 
the kinetic energy. A regular, continuous pattern 
of nodes and antinodes in the kinetic energy exists 
both at the resonant modes and between them. 

(Fig 4). Thus the non-adiabatic pressure work is 
due primarily to the convection. 
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Fig. 4. — Non-adiabatic pressure spectrum at the 
a depth of 100 km. Note that it is featureless 
even at the frequencies of the simulation modes. 
Hence, the non-adiabatic pressure is primarily due 
to random convective processes. 

The mode compression is calculated from the 
Christensen-Dalsgaard modes (since these are es- 
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Fig. 5. — Mode mass in the simulation domain 
and the Sun. Mode mass decreases with increas- 
ing frequency because higher frequency modes are 
more concentrated toward the surface than low 
frequency modes. Because the solar modes extend 
below the bottom of the simulation domain, they 
have a larger mode mass. Because the simula- 
tion domain is shallow, the mode mass becomes 
nearly constant at low frequency where the eigen- 
functions become nearly constant within the depth 
of the simulation. 

sentially identical with the simulation modes at 
the resonant frequencies) normalized by the square 
root of their energy in the simulation domain. Be- 
cause the simulation domain is shallow, while, es- 
pecially at low frequency, the solar modes have 
substantial amplitude below the computational 
domain, the mode mass in the computational box 
is significantly smaller than the actual mode mass 
of the solar modes (Fig 5). 

Finally, we compare the actual mode energy de- 
cay rate with the predicted rate of work by convec- 
tion on the modes given by eqn (5) in figure 6. The 
squares are the actual mode energy decay rates for 
the three resonant modes of the box. The solid line 
is the running mean of the predicted work and the 
dashed line is a smooth two power law fit to the 
predicted work. The agreement is very good, but 
not perfect. There are significant deviations from 
the power law fit in the neighborhood of the modes 
and these are reflected in the actual mode decay 
rate. This is a phenomena that still remains to be 
explained. Notice also that the decrease in work 
toward low frequencies is much slower than for the 
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Fig. 6. — Rate of stochastic energy input to the 
simulation modes (squares) compared to the pre- 
dicted excitation rate (small crosses) from the 
work of non-adiabatic pressure fluctuations on the 
modes (eqn 5). The solid line is a running box- 
car average and the dashed line a two power law 
fit over the entire frequency range. The excita- 
tion rate is larger than the solar case because the 
mass of the modes in the simulation is less than 
the mass of the modes in the Sun. The near con- 
stancy of simulation mode mass at low frequency 
leads to a much slower decrease of excitation with 
decreasing frequency than occurs in the Sun. 

Sun (Fig 7). The reason is the near constancy of 
the mode mass within the simulation domain at 
these low frequencies compared with the steeply 
rising mode mass with decreasing frequency on 
the Sun. This application of our excitation rate 
formula to the modes that are excited in our sim- 
ulation verifies that the formula is correct and can 
be applied to the Sun. 

5. Mode Excitation: Sun 
5.1. Excitation Spectrum 

The excitation rate for solar modes as a func- 
tion of frequency is shown in Fig. 7. To obtain 
these results we used the shorter but higher res- 
olution 125 2 x 82 simulation because it has more 
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Fig. 7. — Rate of stochastic energy input to modes 
for the entire solar surface (simulation=triangles, 
observations=squares) from Roca Cortes (1999), 
based on observed mode velocity amplitudes and 
line widths from GOLF, for modes with t — — 3, 
which are all nearly radial close to the surface. The 
rate of energy input to the solar modes is smaller 
than to the simulation modes by the ratio of the 
mode mass of the solar modes to the mode mass 
in the simulation domain (fig 5). 

high frequency power. The magnitude and fre- 
quency dependence we find for our 6 Mm square 
box is very close to the observed values for the 
entire sun. This means that the pressure fluctu- 
ations must be uncorrelated on larger scales, so 
there is no extra driving contribution. 

What produces this frequency dependence? 
The low frequency behavior is controlled by the 
nature of the eigenmodes and the high frequency 
behavior is controlled by the nature of convection. 
The work integral (Eqn. 5) contains a factor per- 
taining to the modes: the radial gradient of the 
mode displacement normalized by the square root 
of the mode energy (which makes it independent of 
the mode amplitude) . This is small at low frequen- 
cies and increases with frequency approximately as 
^3.5 ^pjg p ar ^ £ ^- ms dependence is due to the 
radial gradient of the displacement. The mode 
dispersion relation is co = y/(gk), so k = u> 2 /g, 
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Fig. 8. — The mode factor in the work integral: 
^f/E^ 2 . Excitation decreases at low frequency 
due to mode behavior, in part because the radial 
wave vector is approximately k — u) 2 /g and in part 
because the mode mass increases with decreasing 
frequency (Fig 5). 

which accounts for two powers of the frequency. 
The remainder of the frequency dependence is due 
to the mode mass, which decreases with increas- 
ing frequency because higher frequency modes are 
more concentrated toward the surface (Fig. 5). 

The mode excitation decreases at high fre- 
quency because the pressure fluctuation power 
decreases with increasing frequency roughly as 
v~ A (Fig. 9). Convective motions, whose power 
decreases at small scales and high frequencies, 
produce the gas and turbulent pressure fluctua- 
tions. This then leads to a similar high frequency 
decline in the pressure spectrum. 

Total pressure fluctuations are small at low fre- 
quency because the atmosphere is nearly in hydro- 
static equilibrium. However, the turbulent pres- 
sure fluctuations are largest at low frequency be- 
cause they are a convective effect and convection 
is a longer time scale phenomena. Hence, the gas 
pressure fluctuations must also be large and out of 
phase with the turbulent pressure fluctuations at 
low frequency in order to produce small total pres- 
sure fluctuations (Fig. 10). The turbulent pressure 
is small compared to the gas pressure (~ 15%), 
but the magnitude of fluctuations are comparable 
for the gas and turbulent pressures, so they can 
indeed cancel each other. 



Z = 200. km 
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Fig. 9. — Spectrum of pressure fluctuations at 
a depth of 200 km, smoothed with running box- 
car mean. The non-adiabatic gas pressure fluc- 
tuations exceed the non-adiabatic turbulent pres- 
sure fluctuations by about a factor of four, but 
they become comparable in the peak driving range 
at larger depths. The local maxima in the to- 
tal pressure fluctuations at 2.6, 3.9 and 5.5 mHz 
are due to the resonant modes in the simulation. 
The non-adiabatic pressure varies smoothly across 
these resonant frequencies, indicating it is primar- 
ily due to convection. The pressure fluctuation 
power decreases roughly as v~ A at high frequency, 
because it is due due to stochastic convective mo- 
tions which decrease in power at high frequency 



5.2. Excitation Location 

At what depth does the driving occur? Con- 
sider the integrand of the work integral at differ- 
ent frequencies (Figs. 11 and 12). At low frequen- 
cies the integrand amplitude is similar over an ex- 
tended depth range, but it is small. Where the 
integrand is large, in the region of peak driving 
around 3-4 mHz, the integrand becomes concen- 
trated very close to the surface and most driv- 
ing occurs between the surface and 500 km depth. 
At still higher frequencies the integrand becomes 
small again and even more concentrated near the 
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Fig. 10. — Correlation of turbulent and gas pres- 
sure at the surface. The turbulent pressure is only 
about 1/6 the magnitude of the gas pressure near 
the surface but the magnitude of their fluctuations 
are similar. 
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11. — Integrand of the work integral, 



w dz 



/8AvE w , (eqn 5) at 2, 3, 4, and 5 
mHz as a function of depth. At frequencies where 
the driving is large, the integrand is significant 
only within 500 km of the surface. 

surface. 

Where in space does this driving occur? The 
warm granules have only small non-adiabatic pres- 
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Fig. 12.- 
grand, u> 2 



Logarithm (base 10) of the work inte- 
nt 2 

SP*§f /8AvE u , (eqn 5) (in units of 

ergs/cm 2 /s) as a function of depth and freqency. 
The work is concentrated close to the surface in 
the peak excitation range (3-4 mHz) and at 
higher frequencies. 




Fig. 13. — Correlation of non-adiabatic pressure 
with vertical velocity at the surface. The large 
negative fluctuations in the non-adiabatic pressure 
occur where the velocity is positive (downward) in 
the intergranular lanes. 



sure fluctuations. Large, negative fluctuations are 
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Fig. 14. — Image of the non-adiabatic total pres- 
sure fluctuations at depth 100 km, with the con- 
tours of zero surface velocity to outline the gran- 
ules. The units of the pressure are 10 3 dynes/cm 2 . 
Large, negative pressure fluctuations occur in the 
intergranular lanes. Positive pressure fluctuations 
are half as large and occur inside the granules. 
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Fig. 15. — Image of non-adiabatic pressure fluctu- 
ations at 100 km depth in the 2-5 mHz range, 
with the contours of zero velocity at the surface 
to outline the granules. The units of the pressure 
fluctuations are 10 3 dyne cm~ 2 . In this frequency 
range, where the driving is maximal, the largest 
SPnad occur at the edges of granules and inside 
the intergranular lanes. 
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Fig. 16. — Rate of stochastic energy input to 
modes for the entire solar surface, showing the 
individual contributions of the non-adiabatic gas 
and turbulent pressure to the work of the total 
non-adiabatic pressure. Most of the driving in the 
peak driving range comes from the turbulent pres- 
sure. 

concentrated in the downdrafts (Figs. 13, 14). 

The maximum mode driving occurs in the fre- 
quency range of 3 - 4 mHz, by non-adiabatic pres- 
sure fluctuations in the same frequency range. By 
filtering the time sequence of these fluctuations 
we see (Fig 15) that in the peak driving range 
also the driving occurs predominantly in the in- 
ter granule lanes and near the edges of granules. 
The high frequency power near granule edges is 
due to the motion of the granule boundaries over 
the one hour time interval on which the filtering 
was performed. This is, in part, a result of changes 
in granule size as they evolve. No direct correla- 
tion of non-adiabatic pressure fluctuations in the 
range of 2-5 mHz with velocity is seen. Keep in 
mind, however, that a correlation plot does not 
reveal correlations of events that happen in the 
neighborhood of one another. 
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Fig. 17. — Mode driving at 4 mHz evaluated from 
the surface to depth z, showing the individual con- 
tributions of the non-adiabatic gas and turbulent 
pressure. Close to the surface the contributions of 
the two are comparable, but there is little contri- 
bution from the gas pressure below 200 km depth, 
while the turbulent pressure work is significant 
down to 500 km depth. 
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Fig. 18. — Horizontally averaged non-adiabatic 
pressure at the surface and emergent intensity 
variation in time. They arc tightly correlated indi- 
cating that radiative cooling at the surface is the 
source of the non-adiabatic pressure fluctuations 
there.. 

5.3. Excitation Source 

What is the source of the non-adiabatic pres- 
sure fluctuations? Is it entropy fluctuations or 
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Fig. 19. — Divergence of the convective and ra- 
diative fluxes, at z=100 km, multiplied by — 
1) compared to the time derivative of the non- 
adiabatic pressure at the surface. The units are 
10 3 ergs/cm 3 /s. The rate of change of non- 
adiabatic pressure closely follows the divergence 
of the net flux but has a slightly larger amplitude. 



Reynolds stresses? Both play a role, but the pri- 
mary source of mode driving is turbulent pressure 
fluctuations (Reynolds stresses) (Fig 16). This 
is surprising, since the non-adiabatic gas pressure 
power is larger than the turbulent pressure power 
near the surface (Fig 9). The non-adiabatic gas 
and turbulent pressures contribute comparably to 
the work near the surface, but the contribution 
of the turbulent pressure extends deeper and pro- 
vides the dominant contribution to the total work 
(Fig 17). 

The gas pressure fluctuations arc significantly 
larger than the turbulent pressure and have max- 
ima at the frequencies of the three radial modes of 
the simulation. The non-adiabatic pressure fluctu- 
ation power, however, varies smoothly across these 
frequencies indicating that it is primarily due to 
stochastic convective processes (Fig 9). 

There is a tight correlation between the non- 
adiabatic pressure fluctuations at the surface and 
the emergent intensity (Fig 18), which indicates 
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Fig. 20. — Correlation of the divergence of the 
convective and radiative fluxes, at z=100 km, mul- 
tiplied by — 1) with the time derivative of 
the non-adiabatic pressure at the surface. The 
units are 10 3 ergs/cm 3 /s. The close correlation 
shows that the slight instantaneous imbalance in 
the radiative and convective flux divergences is the 
primary source of non-adiabatic pressure fluctua- 
tions. 

that it is the fluctuating cooling at the surface 
that is the main source of stochastic mode exci- 
tation there. Indeed, the source of entropy fluctu- 
ations is the cooling of fluid that approaches op- 
tical depth unity (Stein & Nordlund 1998). This 
correlated noise is believed to be responsible for 
the difference in asymmetry of the modal power 
spectra observed in velocity and intensity (Nigam 
et al. 1998; Kumar & Basu 1999). Our discussion 
in terms of non-adiabatic pressure fluctuations is 
equivalent to the discussion in terms of entropy 
fluctuations by Goldreich et al. (1994). 

We can use the energy equation to determine 
the processes producing the non-adiabatic gas 
pressure fluctuations (Stein & Nordlund 1991). 

r P / dP\ 
5P naA = 5P--y5p= [ — J 5s = (T 3 -l)SQ, 

(10) 




Fig. 21. — Logarithm of the gas pressure fluctu- 
ations scaled by the square root of the density, 
as a function of depth and frequency. The color 
scale is identical to the following image of the tur- 
bulent pressure fluctuations. With increasing fre- 
quency the gas pressure fluctuations become more 
rapidly concentrated near the surface and decrease 
more rapidly in strength than the turbulent pres- 
sure fluctuations. 
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Fig. 22. — Logarithm of the turbulent pressure 
fluctuations scaled by the square root of the den- 
sity, as a function of depth and frequency. The 
color scale is identical to the preceding image of 
the gas pressure fluctuations. The turbulent pres- 
sure fluctuations extend deeper and decrease less 
rapidly in magnitude with increasing frequency 
than the gas pressure fluctuations. As a result, 
turbulent pressure is the primary source of p-modc 
excitation in the peak driving range of 3-4 mHz. 
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Fig. 23. — Integrand of the stochastic work inte- 
gral for non-adiabatic gas pressure only, at 4 mHz, 
as a function of depth and time (right panel, eqn 5) 
and the two terms that contribute to it: the hori- 
zontally averaged, non-adiabatic gas pressure fluc- 
tuations, multiplied by (p) -1 / 2 (left panel) and the 
coherent mode density fluctuations multiplied by 
(p) 1 / 2 (center panel). The non-adiabatic gas pres- 
sure fluctuations (left panel) are produced by the 
imbalance of the convective and radiative energy 
transport and are largest close to the surface. The 
color scale goes from maximum negative to maxi- 
mum positive value for each variable. 

and 

^ = ~ (F iad + F conv ) + u-VP + Q diss , (11) 
Dt oz 

where u = u — (pu)/(p) is the convective com- 
ponent of the velocity. The divergence of the 
radiative and convective fluxes are the dominant 
terms, but they nearly cancel each other since en- 
ergy transport shifts from convective to radiative 
near the surface. It is their slight instantaneous 
imbalance locally that leads to the non-adiabatic 
pressure fluctuations (Figs 19 and 20). 

The separate contributions of the gas and tur- 
bulent pressure to the total non-adiabatic pressure 
spectrum are shown in figs. 21 and 22, with the 
same color scale in both. Near the surface, there 
is clearly more power in the gas pressure. How- 
ever, in the peak driving range of 3-4 mHz, there 
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Fig. 24. — The same as figure 23, but for the non- 
adiabatic turbulent pressure only. The turbulent 
pressure fluctuations are more coherent in depth 
and time than the non-adiabatic gas pressure fluc- 
tuations so they produce a more coherent contri- 
bution to the work integral. 

is more power in the non-adiabatic turbulent pres- 
sure below the surface. 

Another reason for the dominance of turbu- 
lent pressure in the mode driving is revealed in 
figs. 23 and 24. The left panels in each show the 
non-adiabatic gas and turbulent pressures (scaled 
by p -1 / 2 ) as a function of time and depth. The 
turbulent pressure varies more slowly with depth 
and has a longer time scale than the gas pres- 
sure. The middle panels show one realization 
of the mode density or compression (scaled by 
p 1 / 2 ). The right panels show the work integrand, 

1/2 

SPned^/dzE' 1 '. The gas pressure contribu- 
tion to the integrand is more concentrated at the 
surface and alternates sign with time and depth 
more rapidly than the turbulent pressure contri- 
bution. The greater coherence of the turbulent 
pressure allows it to make a greater contribution 
to the net work. 

The very largest pressure perturbations are as- 
sociated with the sudden initiation of downdrafts 
and produce waves that propagate up into the 
chromosphere and steepen into shocks (Skartlien 
1998; Skartlien et al. 2000). These events may cor- 
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Fig. 25. — Integrand of the stochastic work inte- 
gral as a function of depth and time (right panel, 
eqn 5) and the two terms that contribute to it: 
the horizontally averaged, non-adiabatic pressure 
fluctuations, multiplied by (p) -1 / 2 (left panel) and 
the coherent mode density fluctuations multiplied 
by (p) 1 ^ 2 (center panel). The non-adiabatic pres- 
sure fluctuations (left panel) are produced by the 
imbalance of the convective and radiative energy 
transport and are largest close to the surface. The 
color scale goes from maximum negative to maxi- 
mum positive value for each variable. 

respond to the large individual acoustic events ob- 
served by Rimmele et al. (1995) and Goode et al. 
(1998). They may be the tail of the distribution 
of the stochastic driving process. A comparison of 
observed acoustic events and those in the simula- 
tion has been made by Goode et al. (1998). 

6. Summary 

How does this all fit together? Figs 25 - 27 
show the time variation of the horizontally aver- 
aged non-adiabatic pressure fluctuations, the den- 
sity (compression) in the modes and the product 
of the two which is the integrand of the work in- 
tegral (eqn 5) as a function of depth and time for 
modes with frequencies 1,3, and 5 mHz. 

The dominant contribution to the work comes 
from the turbulent pressure (Reynolds Stress). 
Near the surface non-adiabatic gas pressure fluc- 



Fig. 26. — Same as Fig 25, but for v = 3 mHz. At 
higher frequencies the mode compression (center 
panel) and hence the work (right panel) becomes 
more concentrated toward the surface. 
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Fig. 27. — Same as Fig 25, but for v = 5 mHz. The 
compression and work are even more concentrated 
toward the surface and the mode has several radial 
nodes. 

tuations, produced by an instantaneous imbalance 
between the divergences of the radiative and con- 
vective fluxes, also contribute. Their divergence 
is individually large, but since energy transport 
switches between convection and radiation in the 
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surface layers, they are nearly equal and opposite. 
At each instant they do not exactly cancel, which 
leads to heating and cooling and hence entropy 
fluctuations. Excitation is small at low frequen- 
cies due to mode properties - the mode compres- 
sion decreases and the mode mass increases at low 
frequency. In addition, at low frequency the mode 
amplitude is small near the surface where the non- 
adiabatic pressure fluctuations are large. (At very 
low frequencies, driving occurs deeper than the 
simulation domain.) Excitation is small at high 
frequencies due to the pressure fluctuation spec- 
trum - pressure fluctuations become small at high 
frequencies because they are due to convective mo- 
tions which have a longer time scale than the dom- 
inant p-mode periods. Large non-adiabatic pres- 
sure fluctuations occur primarily in the intcrgran- 
ular lanes (due to turbulence and surface radia- 
tive cooling) and near the edges of granules (due 
to granule expansion and collapse). 
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A. Evaluation of the Mode Excitation Rate (eqn 5) 

The mode excitation rate (eqn 5) can be evaluated in two different ways: in Fourier space and in real 
space. We describe both methods. The variables that appear in eqn 5 for the mode excitation are the non- 
adiabatic total pressure fluctuation, the derivative of the mode displacement and the mode energy. These 
quantities are calculated as follows for both methods of evaluating the excitation. 

The non-adiabatic pressure fluctuations are obtained from the simulation results. First the gas pressure, 
turbulent pressure (Pt = (pu 2 z )) and density are averaged over horizontal planes and saved at 10 s or 30 
s intervals. These are then interpolated to the Lagrangian frame, that moves with the vertical (radial 
oscillation) motions. The Lagrangian frame is determined by calculating the mass column density at each 
time and interpolating the variables to the time average mass column density. Next, the non-adiabatic total 
(gas plus turbulent) pressure is calculated from 

pnad = (ln( p gas + p tuih) _ Ti lnf)) (p^ + p turb) _ (A1) 

Finally, the fluctuation of the non-adiabatic total pressure about its time average, 

SP nad (z, t) = P nad (z, t) - (P nad (z, t)) t , (A2) 

is determined. 

The mode displacement £ v (z) for the radial mode of frequency v is obtained from the cigenmode calcula- 
tions of Christensen-Dalsgaard, using his spherically symmetric model S (Christensen-Dalsgaard et al. 1996). 
Alternatively, one could extract modes directly from the numerical simulations, using Fourier decomposition. 
A drawback with this method is, apart from slightly noisier mode structure, that there are only three reso- 
nant modes in the shallow simulation box — the 35 radial modes from standard solar models provide much 
better frequency coverage. As we have shown (fig 2), the modes from the simulation are essentially the same 
as the solar modes at the resonant frequencies of the computational domain, where the solar modes have a 
node at the depth of the bottom of the computational domain. 

The mode energy is calculated from the displacement according to eqn. 6. The mode displacements are 
interpolated to the simulation grid at the same total pressure as in the Christensen-Dalsgaard model and the 
derivative of the displacement is calculated. Since the derivative of the mode displacement always appears 
normalized by the square root of the mode energy, the mode amplitude normalization cancels. 

In Fourier space, the work integral which appears in the energy input rate to the modes (eqn 5) is 
evaluated by taking the time Fourier transform of the total non-adiabatic pressure fluctuations at each 
depth, SP nad (z,t). For each frequency and depth, this is multiplied by the spatial derivative of the mode 
displacement at that depth interpolated to the frequencies of the Fourier transform (and normalized by the 

— 1/2 

square root of the mode energy), w(<9£ w / dr)E u . The spatial dependence of the modes varies slowly and 
continuously with frequency, so such interpolation is possible. This product is integrated over the depth of 
the simulation domain for each frequency. The energy input rate to the modes is the square of the absolute 
value of this work integral, divided by the frequency interval for the Fourier transform, Ai>, (which equals 
multiplying by the time interval of the simulation), multiplied by the area of the simulation (36 Mm 2 ) and 
divided by eight (Eqn 5). 

In real space, the integrand of the work integral for each frequency, u, is calculated as the non-adiabatic 
total pressure fluctuation at each depth and time, <5P nad (z, t), multiplied by the normalized derivative of the 

1/2 

mode displacement for that frequency at each depth, uj(d^(z) / dz) / Ej , multiplied by a phase function 
cos (4> + cot) for each time, t, a snapshot was saved in the simulation (10 s or 30 s intervals). For each depth, 
this is summed over all saved snapshots in the longest interval that is an integral number of mode periods for 
the given mode frequency v, and divided by the number of snapshots summed over to get the average value. 
This integrand is then integrated over depth. This is done for two values of the phase, <f> = and tt/2 (since 
the phases between the modes and the pressure fluctuations are random, we average these two orthogonal 
cases). The value of the integral for each phase is squared and these two values are summed. The energy 
input rate to the mode at this frequency is this average multiplied by the time interval integrated over (for 
that frequency), multiplied by the area of the simulation (36 Mm 2 ), and divided by eight. 
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